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Abstract 

We present a unified mathematical approach to epidemiological models with parametric hetero- 
geneity, i.e., to the models that describe individuals in the population as having specific parameter 
(trait) values that vary from one individuals to another. This is a natural framework to model, e.g., 
heterogeneity in susceptibility or infectivity of individuals. We review, along with the necessary 
theory, the results obtained using the discussed approach. In particular, we formulate and ana- 
lyze an SIR model with distributed susceptibility and infectivity, showing that the epidemiological 
models for closed populations are well suited to the suggested framework. A number of known 
results from the literature is derived, including the final epidemic size equation for an SIR model 
with distributed susceptibility. It is proved that the bottom up approach of the theory of hetero- 
geneous populations with parametric heterogeneity allows to infer the population level description, 
which was previously used without a firm mechanistic basis; in particular, the power law trans- 
mission function is shown to be a consequence of the initial gamma distributed susceptibility and 
infectivity. We discuss how the general theory can be applied to the modeling goals to include the 
heterogeneous contact population structure and provide analysis of an SI model with heterogeneous 
contacts. We conclude with a number of open questions and promising directions, where the theory 
of heterogeneous populations can lead to important simplifications and generalizations. 

Keywords: SIR model, heterogeneous populations, distributed susceptibility, final epidemic size, 
heterogeneous contact structure, power law transmission function 
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1 Introduction 

The real-world populations are heterogeneous. The populations consist of individuals, and all the in- 
dividuals are different. This is a basic fact, which does not require any proof. Individuals can differ in 
their age, spatial location, social habits, genome compositions, etc. Inasmuch as we aim to model the 
dynamics of interacting populations, it is an obvious condition to include this heterogeneous structure 
into the mathematical models. There are different kinds of population heterogeneity; here we only 
inspect heterogeneity in population parameters (such as, e.g., susceptibility to a specific disease); the 
parameters are considered as an inherent and invariant property of individuals, whereas the parameter 
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values can vary between individuals. We call such heterogeneity parametric heterogeneity ([23]). Tak- 
ing into account the parametric heterogeneity yields important changes for the population dynamics. 
Most importantly, it means that both evolutionary and ecological aspects of the dynamics have to 
be accounted for. In other words, we must follow not only the total population numbers (ecological 
aspects) but also the changes of the population structure (evolutionary aspects), which are described 
by the transformation of the parameter distribution in the population with time. 

The simplest way to include the parametric heterogeneity into mathematical models is to divide 
populations into subgroups, such that each subgroup has its own specific parameter value. By increas- 
ing the number of groups we can assume that there is a continuous distribution of the parameter in 
the population. 

To illustrate the mathematical approach to the parametric heterogeneity and present the notations 
used throughout the text, we start with the simplest possible mathematical model: the Malthus 
equation 

N = mN, (1.1) 

where ./V = N(t) is the total population size at time t. Equation (jl.ip has the solution N(t) = N(0)e mt 
for the given initial condition iV(0). 

Equation (|l.lj) has parameter m, which is the per capita growth rate (and is often called the 
Malthusian parameter). Model (jl.ip can be used to describe, e.g., the growth of a bacterial colony. 
Employing equation (11. ip as a mathematical model to describe the changes of the population size 
notwithstanding, it is customary assumed that this parameter m is the same for any individual in 
the population. This can be hardly true for any realistic population. The simplest way to describe 
the population which has different values of m is to suppose that the total population consists of k 
subpopulations, each of which has its own parameter value; i.e., we replace (jl.ip with the following 
system of equations: 

Ni = miNi, i = l,...,k, (1.2) 

where now the total population size is N(t) = ^j-/Vj(i). There are k initial conditions for (jl.2p : 
iVj(O), i = 1, . . . , k. Is it possible to describe the behavior of the total population size, whose evolution 
is governed by (|1.2p . using only one equation? Summing all the equations in (jl.2p . we obtain 

N = Y^miNi = N^2 mi ^ = E t [m]N = m(t)N, (1.3) 

i i 

where a natural notation used for the mean parameter value in the population at the time moment t: 
fh(t) = E t [m]=Y / mm(t), Pi (t) = ^, £>(*) = 1. 

Equation (|1,3|) is very similar to (jl.ip . the difference is that it depends on the current parameter 
distribution in the population. If it is possible to find E t [m], then the general solution to (jl.3p is 
simply 

N(t) = N(0) exp{ f E T [m] dr} (1.4) 
Jo 

for the initial population size N(0). 

An important disadvantage of the subgroup approach is that the heterogeneity within each group 
is not taken into account if a fixed number of groups is considered. As it was discussed, in case k — > oo, 
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we can conjecture that there is a limiting continuous parameter distribution. For the Malthus equation 
(jl.ip this means that we suppose that each individual is characterized by its own parameter (trait) 
value u), such that the density of individuals with the given parameter value is n(t, u), and the total 
population size is N(t) = f n n(t, u) du, where Q is the set of the trait values. Now (|1.3p is replaced 
with 



where dt = J^. For (I1.5P we need the initial condition n(0, oS). Using the notation 



for the probability density function (pdf) of the current parameter distribution and the usual notation 
for the mean parameter value 



we find that the solution to (jl.5p is again given by (jl.4p . If one can calculate E t [m] then the solution 
to (|1.5p is found. 

An important remark is worth spelling out. All the models we consider are deterministic. We 
use the probability theory language to describe evolution of the parameter distributions, however no 
stochastic effects are included into our models. 

Models of the form (II. 5p are infinite dimensional systems that describe the evolution of the param- 
eter distribution with time along with the total population size. Such (and more complex) models were 
treated in [3J [5] from the general point of the theory of differential equations in infinite dimensional 
spaces. However, it turns out that many such models can be reduced to low dimensional systems 
of ordinary differential equations (ODEs) [37]. The theory of heterogeneous populations, outlined 
in [37], provides the conditions when the mean parameter value Et[-] can be effectively calculated at 
any time moment using only the knowledge of the initial parameter distribution p(0,u). Therefore 
the heterogeneous models, of which (jl.5p is the simplest example, can be studies analytically. The 
examples include the analysis of heterogeneous Malthus equation (jl.5p [35l 136] , Lotka-Volterra system 
[50J (see also [4j ) , tumor cell dynamics [3D] , and the replicator equation [39] . 

A very rich field for the mathematical models with parametric heterogeneity is theoretical epi- 
demiology, where it is natural to assume that individuals vary with respect to their susceptibility 
to a disease, infectivity to pass a pathogen, contact number, etc. (see, e.g., [6j HU [151 EH EE EH 
[24| [25l l26l l49l I60|). For many models from the cited literature, the general approach of the theory 
of heterogeneous populations yields important analytical results [511 [52] . It is the main goal of the 
present manuscript to review and generalize these results. 

The rest of the paper is organized as follows. In Section [2] we present the necessary mathematical 
facts from the theory of heterogeneous populations with parametric heterogeneity. Section [3] is de- 
voted to formulation of various mathematical models of the epidemic spread that include parametric 
heterogeneity. In Section |4j the machinery from Section [2] is applied to the models from Section O 
Finally, Section [5] is devoted to discussion of some open problems and conclusions. 



dtn(t,uj) = m(uj)n(t,uj) 



(1.5) 
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2 Mathematical theory of heterogeneous populations 



In this section we present the mathematical development of the theory of heterogeneous populations 
with parametric heterogeneity as it is required for the subsequent application to the mathematical 
models of the epidemic spread. We aim to present neither the most abstract possible formulations 
for such models (dubbed systems with inheritance, see [28, 29J) nor the most general form of the 
dynamical systems that can be analyzed by similar means [37] . We refer the reader to the literature 
for a notably wider account of the necessary theory \28 \ \37\ [38] . and remark that our presentation has 
a primary goal to make the exposition self-contained. 

Let us consider two interacting populations, each of which is characterized by its own parameter 
(trait) value at any given time moment, such that we can speak of the parameter distribution in 
each population (the populations are heterogeneous with respect to a given parameter). Denote the 
densities of the populations as m(t,u)i) and n 2 (t,uj2}- This implies that the total population sizes are 
given by 



Ni(t) = / ni(t,uji)duji, N 2 (t) = / n 2 (t,oj 2 ) dui 2 , 

respectively. Here fij, i = 1, 2 are the sets of admissible parameter values, and we suppose that these 
sets are such that the corresponding integrals are always well defined (in fact, in the examples we 
usually suppose that Qj, i = 1,2 are intervals of Mr, but the general theory does not require it). 
An obvious generalization of the considered situation is that k populations are considered, some of 
which could be characterized by their own vector-parameters. Such generalization would require only 
additional notation, and therefore we deal with, without loss of generality, the case of two populations. 
For these two populations the pdfs are given by 

rii(t,u)i) . 
= , 1 = 1 ' 2 ' 

such that pi(t,0Ji) > 0, j n pi(t,uji) duj{ = 1, i = 1, 2 for any t > 0. 

To describe the dynamics of the densities of the populations under study, we start with the general 
model 



d t ni(t,uji) 

m(t,ux) 

d t n 2 (t,uj 2 ) 
n 2 {t,u)x) 



Fi{nx(t,uji),n2(t,u)2)) , 

(2.1) 

F 2 {ni(t,uji),n 2 (t,uj 2 )), 



where F±, F 2 are given functions (the per capita growth rates of the interacting populations) and 
dt = §1 ■ For (|2.ip the initial conditions are ni(0,u;i) and n 2 (0, uj 2 ). 

It is the form of F\ and F 2 that defines the dynamics of interacting populations, and to apply the 
theory of heterogeneous populations, as presented in, e.g., [37], this form must satisfy some additional 
requirements. In particular, we adopt that 

Fi(ni(t,ui),n2(t,uJ2)) = fi{v) + (p(u}i)gi(v), i = 1,2, (2.2) 

where j \, gi, tpi, i = 1,2 are given functions, 

v = (iV 1 ,iV2,^i(i) ) ^ 2 (t)), (2.3) 
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and 

<Pi(t) = E t [ipi] = / (pi(u)i)pi(t,u)i)duji, i = l,2. (2.4) 

In words, we require that the per capita growth rates depend explicitly only on the parameters, total 
population sizes and the mean values of the functions of the distributed parameters. Assume that we 
are interested only in the dynamics of the total population sizes. Integrating the first equation in 
with respect to uj\ and the second one with respect to UJ2, and using (|2.2p . (I2.4p . we obtain 



Ni = Ni(fi(v) + <pi(t)gi(y)), i = 1,2. (2.5) 

Inasmuch as v depends only on the total population sizes and <fi(t), i = 1,2, the dynamics of the 
populations sizes can be found if <fii(t), i = 1, 2 are known. We remark that <Pi(t), i = 1,2 depend on 
the current parameter distributions, which actually has to be found. The remarkable fact is that it is 
possible to calculate (fi(t), i = 1,2 for any time moment if the initial parameter distributions 

are known and two additional differential equations are allowed. This is the main result for our 
presentation. To give it a precise statement, introduce the notations 

Mi(t,X)= I e^*^pi{t,Ui)dui, i = 1,2, (2.6) 
Jn t 

for the moment generating functions (mgfs) of (pi(tOi), i = 1,2 respectively at any time moment. 
The mgf, it it exists, defines uniquely the given probability distribution. In the mgfs for the initial 
distributions we occasionally suppress the dependence on time: 

Mj(A) = Mj(0, A), i = l,2. 

It is a basic fact that if Mj(i, A), i = 1, 2 are known then the mean values of <pi(uji) can be found by 
simple differentiation: 

<Pi(t) = Et[(pi] = d x Mi(t, A)| A=0 , i = 1,2. 
Now we can state the following main 



Theorem 1. Let the dynamics of two interacting populations be described by system (J2JJ) with (|2.2p . 
(|2.3p . ()2.4p . Consider auxiliary variables qi(t), i = 1,2 that satisfy the following differential equations 



<ii = 9i(v), %(0) = 0, z = 1,2, (2.7) 
where gi, i = 1, 2 are as in ()2.2[) . Then 

m u x\ M t (0,A + g t (t)) MjjX + qm 

Proof. From the first equation of (|2.1|) with (|2.2p we have 

m(t,LOi) = ni(0,wi) exp{ / / x (v) dt + <pi{u)\) I g\(\)dt} = 

Jo Jo (29) 

= ni(0,wi) exp{ / /i(v) dt + ipi(ui)qi(t)} . 
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Using 1(23)1 . we find 



ft 



exp{ / fi(v)dt} ni(0,CJi)exp{v9i(wi)gi(t)}da;i = 
Ar\n\ ex Pi / /i( v ) dt } / Pi(0,wi)exp{v?i(a;i)gi(t)}(iu;i 



(2.10) 



exp{ / /i(v)di}M!(0, 9 i(t)). 
o 



1 '* 



iVi(O) 

Using the definition (|2.6f) for mgf, 

Mi(i,A) = / exp{A^i(o;i)}pi(i, wi) dwi = 

= ,tm / exp{A^i(o;i)}ni(i,a;i)da;i = from ([2J1 
iV iW J ft! 

_ &h(y)dt r 

^(0)^1 (t) y ni 

_ /o/i(v)rf* 
A^l(0)iVi(t) 



(2.11) 



exp{(A + 9i(t))v3i(wi)}pi(0,a;i) d^i 
Mi(0,A + gi(*)). 



Finally, applying (|2.10p and (|2.1ip we obtain (|2.8j) . The same calculations are valid for the second 
population. ■ 

Remark 2. Theorem [1] allows to reduce system (|2.ip with (|2.2[) to the four dimensional system of 
ODEs (|2.5p . ()2.7p . where the mean parameter values are 

<Pi(t) = EM = d x In M,(A)| A=9i(t) , i = 1,2. (2.12) 

It is a simple matter (see, e.g., |37| ) to prove that the mean parameter values satisfy the equations 

&(t)=#(v)0?(t), * = 1,2, (2.13) 

where cf(t), i = 1,2 are the current variances of tpifai), i = 1,2. 

Remark 3. Using Theorem [T] we find that the solution for the total population size of the heteroge- 
neous Malthus equation (|1.5|) is given explicitly as 

N(t) = iV(0)M(i), 

where M(i) is the mgf of the initial parameter distribution p(0,u)). This solution was studied in detail 
in [341 [35]. 

Remark 4. Theorem [1] shows that the mgfs at any time moment can be found from the mgfs at the 
initial moment using (|2.8p . This theorem also shows that by reducing the original infinite dimensional 
system (|2.ip with (|2.2|> to the system of ODEs (12. 5j) . (12. 7p we actually do not loose any information 
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because, as far as the functions qi(t), i = 1,2 are known, the densities nj(t,Wi), i = 1,2 at any time 
moment can be found from the total population sizes and the time dependent probability distributions 
Pi(t,uJi), i = 1,2, which can be inferred from the corresponding mgfs. For the following exposition we 
do not require explicit formulas for the evolution of distributions, and refer the reader to |37| . where 
these results can be found. 

Remark 5. Theorem [T] provides an analytical proof of (I2.8p . It is interesting to note that there is a 
simple probabilistic proof of a similar formula [IJ [2] . 

Consider the so-called proportional frailty model [2]: it is assumed that the hazard rate of an 
individual is given as the product of an individual specific quantity Z and a basic rate a(t): 

a(t\Z) = a(t)Z. (2.14) 

Here Z plays the role of the parameter distributed in the population. Given Z the probability of 
surviving up to time t is 

S(t\Z) = e - ZA ^\ A{t)= f a(r)dr. 

Jo 

The population survival function is therefore 

S(t) = P[T >t} = E[e~ ZA ^] = M{-A{t)). (2.15) 

The frailty distribution in the population surviving at time t can be found as follows (here l(T > t) is 
the indicator function): 

M(t,A) = E t [e^] = E[e**|T > t) = = 

E[e xz - ZA ^} _ M(X-A(t)) (2 ' 16) 
E [ e -ZA(i)] - M(-A(t)) ' 

Using q(t) = -A(t) in (I2TT6D we obtain ([21^) . 

Concluding, using the mathematical theory of heterogeneous populations with parametric het- 
erogeneity, outlined in this section, it is possible to model communities of populations when each 
population is characterized by its own parameter value, and the per capita growth rates depend on 
this parameter and the average characteristics of the interacting populations. 



3 Heterogeneous models in epidemiology: Model formulation 

The modern mathematical epidemiology has its roots in now classical work by Kermack and McK- 
endrick |2H I43j. where the total population of the constant size N was subdivided into three categories: 
susceptible individuals S that are prone to infection, infectious individuals I that transmit the disease, 
and removed individuals R that either acquire life-long immunity or die. In the simplest case, assum- 
ing that the period of infection is exponentially distributed with the mean I/7 and the transmission 
process is described by mass-action kinetics |31| |4"T] (this means that the contact rate is proportional 
to the total population size N), we have that the dynamics is described by 

S = -/3SI, 

i = (3SI- 7 I, (3.1) 
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where f3 is the transmission parameter that encompasses the information on the probability of a 
successful contact (i.e., the contact that results in infection) and the contact rate. From (|3.ip it 
follows that S(t) + I(t) + R(t) = N for any t, therefore the third equation is redundant and in the 
following we usually omit it. SIR model (|3.ip includes two parameters that in reality vary from 
individual to individual. Therefore, it is important to take this heterogeneity into account. 



3.1 Heterogeneous SIR model with distributed susceptibility to a disease 

There are a number of studies in the literature which model heterogeneous susceptibility to a disease, 
with either a finite number of different susceptibility classes [131 ISS EZJ [32j [331 EE EZ] or with 
a continuous distribution of susceptibility [HI [161 [Ml ESI [26] . We present a formulation from [51] 
keeping in mind that both discrete and continuous distributions can be accommodated. 

Let us denote s(t,ui) the density of the susceptible individuals having the trait value u. The total 
size of the susceptible population is S[t) = j n s{t,uj) dw. Assuming the law of mass action, we obtain 
that the changes in the susceptible and infectious populations are described by 



d t s(t,uj) = -/3(u)s(t,u)I(t), 

i(t) = I(t) [ p(w)s(t,u)dw-'yl(t). 



(3.2) 



Or, using the notations 

s(t, (jj) 



3.21) can be rewritten as 



d t s(t,u) = -p(u)s(t,u)I(t), 

i(t) = p(t)si--ri(t). 

For (|3.3p the initial conditions are s(0,ui) = sq(uj) = Sop s (0,oj), 1(0) = Jo- 
Model (13. 2jl . (13. 3p can be also obtained from the general epidemic equation pj) 



(3.3) 



dts(t,uj) = s(t,u) / / A(r,u,rj)dts(t — T,rj) drdrj, (3-4) 
Jn Jo 

where A(r, u},rj) is the expected infectivity of an individual that was infected r units ago while having 
trait value r] towards to a susceptible with trait value u. If we assume that A(r,uj,r]) = P(u)f{r) for 
a given function / and set 

I(t) = ~ I I f(r)d t s(t-T,r,)drdr,, 



n Jo 



then, after some algebra, we obtain 

d t s(t,uj) = -p(u)s(t,u)I(t), 



7(t) = /(O)J(t) / P(0j)s(t,u)dl0- I I f'(T)dt8(jt-T,Tl)dTdTI. " ^ 



n Jo 



Letting f(r) = e 7r reduces (|3.5p to (|3.3p (this can be shown to be the only case to end up with an 
ODE system, for this it is necessary and sufficient that / satisfies the equation / = —7/). 
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Remark 6. If in the general epidemic equation f|3.4[) it is assumed that /(r) = x(T — r), where % is 
the Heaviside function then from (|3.5p it follows that 

dts(t,u)=-p(u)s{t,u)I(t), 

i(t) = j3(t)SI - P(t - T)S(t - T)I(t — T). 

Model (|3.6p is the model with distributed susceptibility studied in [231 [251 12S] ■ 



3.2 Heterogeneous SIR model with distributed infectivity 

Let /3(u) be the transmission parameter of an individual with the infectivity value cj, and i(t,co) be 
the density of the infectious hosts with trait value d at the time moment t, I(t) = f n i(t,uj) duj. We 
assume that the susceptible population is homogeneous. To describe the dynamics of the infectious 
population it is necessary to specify what trait value is assigned to a newly infected individual, which 
was infected by an infectious individual with the trait value rj. Denoting ip(u), rj) the pdf that prescribes 
the probability that a newly infected individual is assigned the trait value u) if infected by an individual 
with the trait value rj, we obtain 

S(t) = -P(t)S(t)I(t), 

f (3-7) 

d t i(t,u) = S(t) / ip(u,r})/3(r])i(t,rf)dr) -ii(t,u), 



n 

where now 

i(t, oS) 



P{t) = / f3(u;)pi(t,uj)du}, Pi(t,u) 



The initial conditions are 5(0) = So, i(0,uj) = io(uj) = Iopi(0,u). 

There are several natural choices for if)(u, rj), the simplest of which is rj) = 5(lo — cj'), where 5 
is the delta-function. This option means that the newly infected individual acquires the trait value of 
the person by whom he was infected (this is equivalent to ip(u),r)) = pi(t,oj), i.e., the trait values are 
assigned according to the current distribution of the infectivity). Using ip(uj,r]) = 5(u — uj') in (13.71) 
we obtain 

S(t) = -P(t)S(t)I(t), 
d t i(t,u) = (3{uj)i(t,uj)S(t) -ji(t,aj). 



Model (13. 8p is very similar in form to (13. 3D . However, it should be clear that, in general, model (13. 3D 
is much closer to reality than (|3,8p . which corresponds to the case when several different strains of an 
infection can be passed on. 



3.3 Heterogeneous SIR model with distributed susceptibility and infectivity 

Let us assume now that both the susceptibility and infectivity are distributed in the population 
experiencing the disease; this will generalize models (|3.3[) and (|3.8p . 

Let s(t,ui) and i(t,ui2) be the densities of the susceptible and infectious individuals respectively. 
For the following, simplifying, we assume that the traits of the two subpopulations are independent, 
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i.e., (3(loi,u)2) = Pi((jJi) 132(1^2) ■ The number of susceptible hosts with the trait value u\ infected by 
infectious individuals with the trait value u)% is given by 



/3i(wi)s(t,Wl)^2(w2)i(t,W2)- 

Therefore, the total change in the infectious subpopulation with trait value 0J2 is 



assuming that ijj(u),r]) = 5{oj — rj) (cf. (|3.8p ). A similar expression describes the change in the 
susceptible population. It is worth emphasizing that nothing else except for the standard law of mass 
action is supposed to formulate the terms for the change in susceptible and infectious subpopulations. 
Combining the above assumptions we obtain the following model: 



8 t s(t,UJl) = -Pt(tjOl)s{t,U)i) I (32(uJ2)i(t,LJ2)duJ2 

Jn 2 

= -f3 1 (u; 1 ) S (t,u Jl )(32(t)I(t), 



d t i(t,u! 2 ) = fh{w2)i(t,u)2) / Pi{ui)s{t,u)i)dui-"ji{t,u)2) 
= #2(w2)i(t, U2)Pi(t)S(t) - ji{t, wa), 

where 

fait) = E t [Pi] = / ,9i(wi)p a (t,wi)ci»/i, Mt) = £t[fo]= f fafa)Pi{t,u} 2 )du 2 . 

Model (|3.9p is supplemented with the initial conditions s(0,u;i) = Sop s (0,uJi), i(0,u>2) = IoPi(0,u>2)- 

In the case 7 = we obtain the heterogeneous SI model with distributed susceptibility and infec- 
tivity 

d t i(t,L0 2 ) = /h(W2)i(t,U} 2 )Pl(t)S(t), 

for which the global dynamics is simple and is similar to the simplest homogeneous SI model, when 
S(t) -4 0, I(t) -4 N, when t -)■ 00. 



3.4 Heterogeneous SI model with heterogeneous contact structure 

Above we discussed mainly about heterogeneity of the hosts: whether all susceptible individuals are 
of the same type with equal susceptibility, and whether all infectious individuals have equal ability to 
infect others. Another aspect of heterogeneity is the possible heterogeneous social contact network, 
which is the one of the central topics in mathematical epidemiology, e.g, \TT\ [T7l H2l I48j . It is 
difficult to apply the general theory of heterogeneous populations as presented in Section [2] to such 
models, however, there is a simple case, for which some results can be obtained. 

Let us assume that n(t,w) denotes the density of individuals in the population, which are making 
ui contacts on average. Every individual can be contacted by another individual, and the individuals 
differ in an average number of contacts. This situation is usually termed as separable mixing. If we 
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denote r the probability of transmission the disease given a contact, then, the simplest SIR-model 
with separable mixing can be described by the following system: 

„ , , f n uji(t,uj)du 

d t s(t,w) = -rus(t,u^ 



f a u}8(t,u)du + f n ui(t,u)dw ' 

dti(t, w) = -rus{t, u)- — — — — — — — ji{t, u) . 

j n ujs{t,u) du + Jq ui{t,uj)du} 

In case of SI model (|3.10p (7 = 0) we have s(t, u) + i(t,uj) = no(u) for any t, and no(u) is a given 
density which specifies probability density function of the contact distribution. Using the property 
i(t,uj) = no(w) — s(t,oj), we obtain 



dts(t,w) = —rujs(t,uj) 



_ f n us{t,u)du 
ojno{uj)duj 



(3.12) 



We conclude this section with an obvious statement that the list of possible epidemiological models 
with parametric heterogeneity can be made longer. What is important, however, that models (|3.3|) . 
(I3~6j) . (I3~8l) . ([531) . (I3TTUD . and (i3~T2l all are written in the form ipi)H [23 )l . which allows a unified 
treatment of these models within the mathematical framework outlined in Section [2J 



4 Mathematical analysis of an SIR model 

4.1 Analysis of heterogeneous SIR model with distributed susceptibility and in- 
fect ivity 

Model (|3.9p comprises, as particular cases, models fj3.3j) and f|3 . 8j> . Therefore, the main result is stated 
for this model (the case of distributed susceptibility alone was treated in [51]). 

Theorem 7. Heterogeneous SIR model (|3.9p lozf/i distributed susceptibility and infectivity is equivalent 
to the following two-dimensional non- autonomous system of ODEs: 

S = -h 1 (S)h 2 (t,I), 

(4.1) 

I = hx{S)h 2 {t,I) -7/, 

where 

/H(S) = SoCSaM^CO, A)!;^!?,,)- 1 , ( 4 - 2 ) 
/ia(t, J) = /oe^(3 A MjH0, AJL^/j,,)- 1 , (4.3) 

and M^O, A), i = 1,2 are ffoe inverse functions to the mgfs of the initial distributions of susceptibility 
and infectivity respectively. 

Proof. According to Theorem [fl system (|3.9p is equivalent to 

s = -h{t)h{t)si, 
i = h(t)h(t)si- 1 i, 

91 = -&(*)!, 
h = M)S, 
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where 



&(i) = d A lnMi(0,A)| A=? . (t)) % 



= 1,2. 



From (I4.4p it follows that 




-lnM 1 (0,q 1 ), 
at 

— lnM 2 (0,g 2 ) -7, 



from which we have two first integrals 



^ = M 1 (0, ?1 ) 



M 2 (0,g 2 ). 



(4.5) 



Since the mgfs are monotone functions, we can express q\ and r/ 2 in (|4.5|) through 5 and / respectively. 
Putting the resulting expressions into (14. 4p and remembering the inverse function theorem, we obtain 



Theorem [7] gives an important example of a bottom up approach in building mathematical models 
in epidemiology [181 130] . We start with a detailed model when each individual in the population has 
its own parameter value. This model, albeit still oversimplified and unrealistic, takes into account 
physiological structure of the population without any ad hoc assumptions on the dynamics of the 
total population sizes. The general theory of the heterogeneous populations with parametric hetero- 
geneity allows to reduce the model to the population-level description, which, incidentally, still can 
be described as a classical SIR model, which contains non-linear and time-dependent transmission 
rates. This shows that a great deal of mathematical models, built from the first principles, are still 
valid and can be used for ecological predictions, because they can be shown to follow from detailed 
individual-based models. 

Let us consider an explicit example (more details are given in |52|). 

On the power law transmission function. The transmission function (the number of new in- 
fectious cases per time unit) is considered to be one of the major ingredients of the models describing 
spread of an infectious disease |19] , Historically, borrowing an analogy from the chemical kinetics, 
this function was supposed to be of a simple bilinear form, oc SI, which means that the law of mass 
action is assumed to hold |31| (this is often called density dependent transmission function). If one 
assumes that the number of contacts is fixed for any individual and does not depend on the popula- 
tion size, than the transmission function takes the form oc SI/N (frequency dependent transmission 
function). If the population size constant these two transmission function yield the same predictions, 
whereas variable population size can produce different behaviors (e.g., |12t I54|). Other transmission 
functions are also possible, see [47] . Important point here is that most of the used nonlinear functions, 
that are somewhat intermediate between density and frequency dependent modes of transmission, are 
phenomenological and lack any mechanistic derivation. 

One of most frequently used transmission function takes the form 



(gZED with (02]), (l4~3l) . 



T(S,I)=P&P, p,q 



> 0. 



(4.6) 
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This function is often called the power law transmission function. It was used in, e.g., |61l I62j in the 
form T(S,I) = j3S p I. Severo [58] considered the full form (|4.6p . see also [IH 05] for a detailed math- 
ematical analysis of epidemiological models with the power law transmission function. Simulations 
|56[ [59] show that power law transmission function indeed improves an accuracy of the mean-field SIR 
models. Using Theorem [7] we can show that the power law transmission function can be inferred from 
the heterogeneous SIR models (|3.9p and (|3,10p . 

Let us assume that the initial distribution of susceptibility is a gamma-distribution with parameters 
k\ and v\\ 

u kl 

p 1 (0,co) = —\—co k ^ 1 w-^, u;>0, ki,u 1 >0. (4.7) 

r(fci) 

From (|4.7p one has, assuming additionally that /3i(ct>i) = uj\, 

M l (\)=(l-±\ \ (4.8) 



therefore the expression in (|4.2p is given by 



h „ ( S 



i/fci 



An analogous expression can be obtained for h,2(t,I) if it is postulated that the initial distribution 
of infectivity is also a gamma-distribution (|4.7p with parameters /c2)^2 and ^2(^2) = Expression 
■9h implies 



Corollary 8. The power law transmission function (14.61) with the heterogeneity parameters q = 1, p = 
1 + 1/ki can be obtained as a consequence of the distributed heterogeneous SIR model (|3.3p with dis- 
tributed susceptibility if the initial distribution of susceptibility is the gamma- distribution with param- 
eters k\, v\. 

The power law transmission function (14. 6p with the heterogeneity parameters q = 1 + 1/&2, P = 1 + 
1/ki can be obtained as a consequence of the distributed heterogeneous SI model (|3.10p with distributed 
susceptibility and infectivity if the initial distribution of susceptibility is the gamma- distribution with 
parameters ki, v\ and the initial distribution of infectivity if the gamma-distribution with parameters 

k2,V2- 

Remark 9. Using (|2.8p and (|2.12p it is straightforward to show that for the initial gamma distribution 
of susceptibility, the time-dependent distribution is also the gamma-distribution with parameters k\ 
and V\ — qi(t), where qi(t) is the solution of the corresponding auxiliary equation (|2.7p . That is, we 
have 

hit) = EM = k ± 1 - , alit) = ^™ . 

v x -q x (t) (^i-gi(t)) 2 

Note that for any time moment the coefficient of variation cv = o~\{t) / f3\{t) = l/^/kl and does not 
depend on t. 

Remark 10. To analytically analyze models with parameter distributions it is useful to have families 
of distributions for which their mgfs are known, as in the case of the gamma-distribution. A very gen- 
eral family of distributions is the so-called power variance function distributions (PVF distributions) 
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[2], which are defined through mgf: 



M(A) = exp 



-P 1 



A 



(4.10) 



with v > 0, m > —1, mp > 0. 

The mgf ()4.10p describes a number of distributions. For instance, if p — > oo and m — > in such a 
way that pm — > k, then (I4.10P reduces to the mgf of the gamma-distribution (14. 8f) . When m > then 
(I4.10p gives a compound Poisson distribution, which is the sum of independent gamma-distributions 
with the parameters v and m, when the number of summand is Poisson distributed with expectation 
p. Such compound Poisson distribution has non-zero mass probability at zero, which means, in terms 
of susceptibility to a disease, that there is part of population of mass exp{— p} that is immune to the 
disease. When m = —1/2 and p < we have an inverse gaussian (Wald) distribution. Other limiting 
distributions are possible, which gives a wide choice of heterogeneity distributions [2]. 

From (I4.10P it can be shown that 

E[Z] = P —, Var[Z] = ^^±i. 

V V V 

Equation (|2.8p implies that the time dependent parameter distribution of the models with parametric 
heterogeneity is still PVF distribution if the initial one is given by (|4.10p . with the following change 
in parameters: 

v-q{t)J 

Not all the distributions possess this "stability" property: e.g., the initial uniform distribution turns 
into truncated exponential distribution. 



On the final epidemic size. A very important quantity of the SIR model (|3.ip is the final epidemic 
size, which can be defined as the number of susceptible hosts that escape infection, and which we denote 
Soq. From f)3. 1 1) it follows that in the homogeneous case this number can be found as the root to the 
equation 

= S e^- N ^\ (4.11) 
Consider again the heterogeneous SIR model (|3.3p with distributed susceptibility. Theorem [7] implies 

Corollary 11. The final epidemic size of the heterogeneous SIR model with distributed susceptibility 
(|3.3p can be found as the root to the equation 



S oo = S M((S oo -N)/ 7 ). (4.12) 
Proof. From the first equation in (14. 4p and recalling the third equation in (|3.ip we find that 

|ln5(t) = |lnM(- J R/ 7 ), 

from which, after integration and using the fact i?oo = N — S^, (|4.12|) follows. ■ 

Remark 12. Equation (|4.1ip is a particular case of (|4.12j) if the initial distribution of susceptibility 
is the delta-function. 
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Figure 1: Final epidemic sizes versus the initial variance of the susceptibility distributions for two 
initial distributions: a gamma distribution (the top curve) and a Wald distribution (the bottom 
curve). All other parameters are the same for both cases 



Remark 13. Equation (14.12P can also be obtained from the general epidemic equation (13. 4j) . which 
shows that no assumptions on the distribution of epidemic length is required to deduce (|4.12p . A very 
careful mathematical analysis of this equation is given in |41| , and we refer the reader to this reference 
for many intricate details. We note that for the first time equation (14. 12f) was written in [9] for a 
discrete distribution of susceptibility. It can be easily proved that the final epidemic size equation 
implies that the epidemic is the most severe in case of a totally uniform population; heterogeneity in 
susceptibility increases the final epidemic size Soo (see Fig. [[]). 

Consider a numerical example. Let us assume that we have two SIR model (|3.3p with heterogeneous 
susceptibility, and the initial distributions are given by a gamma and Wald distributions respectively 
with the same initial mean values and the same variances. The final epidemic sizes are shown in Fig. 
[1] depending on the initial variances. 

Figure [1] allows to make an important conclusion: to infer the final epidemic size for an SIR 
model with distributed susceptibility it is not enough to know only several first moments of the initial 
distribution. This conclusion is somewhat restrictive for a predictive use of such models; however, it 
also signifies that various approximation techniques can lead to erroneous conclusions. As an example 
consider the final epidemic size equation for (|3.6p . obtained in [26] with two different methods. First it 
was supposed that the susceptibility distribution is a gamma-distribution. The second approach was to 
consider an infinite dimensional system of ODE, which can be inferred from (|3.6p . for the moments of 
the corresponding distribution (this is a natural strategy for analyzing infinite dimensional dynamical 
systems of the form (12. 2p . see [491 160]). Inasmuch as it is impossible to solve a system with infinite 
number of equations, various techniques to close such systems exit. In particular, it was supposed 
in [26] that the coefficient of variation is constant. This assumption led to the same result as was 
obtained for the initial gamma distribution. Therefore, it was concluded that the exact form of the 
initial distribution is irrelevant because two seemingly different approaches lead to the same outcome. 
However, theory from Section [2] shows that the opposite is true. First, Remark [9] gives the explanation 
why two approaches in [26] turned out to be equivalent. And second, Fig. Q] provides indirect proof 
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Figure 2: Solutions to the heterogeneous SIR model (|3.3j) with distributed susceptibility. The pa- 
rameters are Sq = 999, Iq = 1, 7 = 700, /?i(0) = 1 and the initial distribution of susceptibility is 
the gamma distribution with erf (0) = 0, 0.2, 0.4, 0.6, 0.8, 1, 1.2 (from the bottom to top curves for the 
susceptible hosts, and the opposite direction for the infectious hosts) 



that the final epidemic size depends on the exact form of the underlying susceptibility distribution 
(see also [41J for exact mathematical statements). 

Unfortunately, for the general case of the SIR model (13. 9p with distributed susceptibility and 
infectivity we were not able to obtain a simple equation for the final epidemic size. 

Numerical illustration. In Fig. [2] numerical solutions for (|3.3p are shown for different initial 
variances of the parameter distribution and equal means, which confirms that the heterogeneity in 
susceptibility decreases the severity an an epidemic. 

In Fig. [3] numerical solutions for system (j3.8H are shown for different initial variances of infectivity 
and equal means. This figure shows that the distributed infectivity has an opposite effect on the sever- 
ity of an epidemic comparing with the distributed susceptibility: the more heterogeneous the initial 
distribution of infectivity is, the smaller the number of susceptible hosts that escape the infection. 

The last conclusion can be cast in an exact mathematical statement if one is interested only in the 
initial period of the epidemic. 

Corollary 14. Let Si(t), S 2 {t) be the solutions to (|3.3ft with the initial conditions that satisfy o~f(0) > 
<r| (0) for the distribution of susceptibility, all other initial conditions being equal. Then there exists 
e > such that Si(t) > S 2 (t) for all t G (0,e). 

Let Si(t), S2(t) be the solutions to f|3.8[) with the initial conditions that satisfy <rf (0) > cr|(0) for 
the distribution of infectivity, all other initial conditions being equal. Then there exists e > such that 
Si(t) < S 2 (t) for allte (0,e). 
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Figure 3: Solutions to the heterogeneous SIR model (|3.8p with distributed infectivity. The parameters 
are So = 999, Iq = 1, 7 = 700, /?i(0) = 1 and the initial distribution of infectivity is the gamma 
distribution with of(0) = 0,0.005,0.01,0.015,0.02,0.025,0.03 (from the top to bottom curves for the 
susceptible hosts, and the opposite direction for the infectious hosts) 



Proof. Using (|2. 13|) and differentiating the first equation in (j3.3[) we obtain 

S"(t) = I 2 S(a 2 {t) + p(t)) - mi'S, 

or, at the initial time moment, S"(0) > S'^O), which proves the first part the corollary. The second 
part is proved in a similar way. ■ 

Finally, for the general heterogeneous SIR model (|3.9|) with distributed infectivity and susceptibility 
some numerical results are shown in Fig. HJ where the initial distributions of susceptibility and 
infectivity are combined from those in Figs. [2]and[3l The conclusion from numerical calculations is 
that the interaction of the distributions of susceptibility and infectivity is nonlinear and cannot be 
predicted from knowing the final outcomes of separate epidemics (see also |8| for the discussion of the 
final epidemic size of the model (13. 9|) ). 



4.2 Heterogeneous SIR model with distributed contact number 

Here we start with system (|3.11|) : 

dts(t,w) = —rus(t,u;)fj,(t), 

dti(t, w) = —rujs(t, oj)fj,(t) — 72 (t, u)) , 

where 

f a ui(t,u)dw 



(4.13) 



Jo ojs(t, u) duj + j n ui(t, ui) cko 
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Figure 4: Solutions to the heterogeneous SIR model (|3.9p with distributed infectivity and susceptibility 
. The parameters are So = 999, Iq = 1, 7 = 700, /?i(0) = 1 and the initial distribution of susceptibility 
is the gamma distribution with of (0) = 0,0.2,0.4,0.6,0.8, 1, 1.2, the initial distribution of infectivity 
is the gamma distribution with a%(0) = 0,0.005,0.01,0.015,0.02,0.025,0.03 



First we note that system (I4.13D is not covered by Theorem [H therefore, at least using the framework 
from Section [21 we cannot reduce this system to an ODE system. However, some results are still 
possible to obtain. In particular, we give a simple derivation of the expression for the initial growth 
of fj,(t) (cf. [35]). 

If we denote K(t) = f n u)i(t,uj), then 



Differentiating fi{t) yields 



K'(t) =rE[u 2 }Sfi{t) -jK(t). 



K'(t) 

At) = 7 77 w .I f — 77 w + * = 

rE[u> 2 ]Su(t) o, . 



(4-14) 



" W '. EH5 + W) " 7 ) +7 " W 
,„ ( )['^- 7 V7f 2 W, 



where the last equality holds for t — >• because S(t) —> N and — )■ 0. In words, we obtained the 
well known results [46] that initially the change in the number of susceptible hosts is proportional not 
to the mean number of contacts but to the mean number plus the coefficient of variation. This shows 



18 



that the individuals who have a high number of contacts contribute disproportionally to the spread 
of an epidemic. 

Corollary [T3] implies that in case of distributed susceptibility or infectivity the initial phase of an 
epidemic can be determined by the first two moments. Here we prove that this is not so for model 
(|4.13|) , as it was shown in [UD] . The proof uses Theorem [H direct proof is given in [BO] . We use the 
notation L(t) = j n ujs(t,u}) du>. 

Lemma 15. Let Si(t), S^i) be the solutions to model (|4.13p with the initial conditions such that 
o~f(t) > o^it) for the distribution of the contact number, all other initial conditions being equal. If 
L(0) < K(0) then there exists e>0 such that Si(t) > S 2 (t) for all t £ (0,e). If L(0) > K(0) then the 
opposite holds. 

Proof. Prom the first equation in (|4.13p it follows 

S" = -rE'[oj]Sfi(t) - /i(t)E[u)]S' - E[u\S/J?(t) = (from and gUD) 



ra 2 (t)n 2 (t)S - fi(t)E[u]S 



a 2 (t)ii 2 (t)S 



1 



Lit) 
K(t) 



rE[u 2 ]S 
L(t) + K{t) 

+ ..., 



+ . . . 



(4.15) 



where dots denote terms that depend only on the first moment of the contact distribution. 

From (I4.15P lemma follows. ■ 

As it was mentioned, Theorem [7] cannot be applied to (|4.13p . However, if we consider the case 
7 = 0, we still can reduce the heterogeneous model to ODE. We rewrite equation ()3. 12|) in the form 



dts(t,w) = —rcjs(t,u) 



1 



<D(t)S(t) 
C 



(4.16) 



where C is the number of contacts, which are made by the total population, u(t) is the average number 
of contacts made by one susceptible individual at time t. Theorem [7] implies that (|4.16p is equivalent 
to the following ODE: 

S(t) = -rh(S) ^ k{S) 



C 



(4.17) 



where h(S) is given by (|4.2|) . For example, if the initial distribution of the number of contacts of 
susceptible hosts is a gamma distribution with parameters k and u, equation (|4.17p takes the form 



S(t) 



s_ 



l/k 



S I 



S_ 

So 




(4.18) 



Numerical solutions of (|4.17p for two parameter sets are given in Figs. [5]and[6j These figures show 
that initially the heterogeneous contact rates increase the speed of an epidemic; in the long term run, 
however, the presence of individuals who make a small number of contacts, slows the epidemic down. 
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Time 

Figure 5: Solutions to the heterogeneous SI model (14.18!) with distributed contact number. The 
parameters are So = 999, 1$ = l,u)(0) = 1 and the initial distribution of contact number is the gamma 
distribution with a 2 (0) = 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 




Time 

Figure 6: Solutions to the heterogeneous SI model (|4.18p with distributed contact number. The 
parameters are Sq = 999, Iq = l,u)(0) = 5 and the initial distribution of contact number is the gamma 
distribution with a 2 (0) = 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 



20 



5 Concluding comments 



In this manuscript we reviewed and generalized several known results on the spread of epidemics in 
a closed heterogeneous population obtained with the help of the theory of heterogeneous populations 
with parametric heterogeneity. Our discussion was mainly about the SIR model (|3,9p and its partic- 
ular cases (|3.3|) . (|3.8|) . and (|3.1(jp . We also introduced the methods of the theory of heterogeneous 
populations to the models that take into account non-uniform contact structure of the population, 
noting that explicit results can be inferred only for a very limited case of the SI model ()4.16p . 

While speaking of heterogeneity of the populations experiencing a disease, three different sources 
of heterogeneity can be accounted for. First, this is the heterogeneity in disease parameters that we 
termed the parametric heterogeneity, and we discussed this at length in the main text. Secondly, 
this is the heterogeneity of the social structure, which results in a different contact rates for different 
individuals. Finally, there is a third source of heterogeneity: the distribution of the infection length. 
For most models in the text it was tacitly assumed that this distribution is exponential with the mean 
I/7 for all the infectious individuals. This assumption is usually made to simplify the mathematics 
and does not correspond to real situations. Using the approach from Section [2] it is possible to 
slightly generalize our models by assuming that parameter 7 is also distributed through the population. 
Proceeding along the lines of Section [3l we obtain 

* = -I" 1 , 1 - (5.1) 
dti(t,ui) = f3Si(t,oj) — 7(o;)i(t, u), 

which is equivalent, according to Theorem [H to the system 

S = -0SI, 

(5.2) 

I = f3SI-d x \nM(\)\ x= _ t I, 

with the mgf M(A) of the initial distribution of 7. Systems (|5.ip and (|5.2p assume that the population 
consists of subpopulations, each of which has an exponentially distributed infection length, but this 
length varies from group to group according to the given initial distribution. A somewhat more 
interesting approach is to consider instead of (|5.ip the following system 

S = -P(t)SI, , . 

(5.3) 

dti(t,uj) = (3(co)Si(t,uj) — ^y(uj)i(t,uj), 

where (3(oj) and 7(0;) are correlated. Model (|5.3p is not covered by Theorem[TJ but still can be tackled 
with the general theory from [37J. 

In the mathematical models considered in the text it was always assumed that the population size 
is closed. This assumption allowed to formulate the models in the form suitable for Theorem [TJ It is 
a natural extension to consider models, when the demography and immigration processes are taken 
into account. This is very important because heterogeneous susceptibility of many diseases can be 
explained by a heritable genetic basis. For example, we can consider the simplest SIR model with 
distributed susceptibility and recruitment in the form 

dts(t,uj) = As(t,u) — P(u))s(t,u)I, 

(5.4) 

I = /3(t)SI — 7 J. 
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Model (|5.4p can be reduced to an equivalent ODE system using Theorem [TJ However, the long-term 
behavior of (|5.4p is straightforward: the individuals with higher values of the susceptibility parameter 
will outcompete those with smaller ones. For the model to be realistic it is also required to include the 
stochastic hereditary element — mutations. That is, a more realistic counterpart of (|5.4p is probably 

dts(t, ijj) = A / ip(uj,r])s(t,rj) drj — (3(u})s(t,u)I, 

Jn (5.5) 

i = P(t)SI- 7 I, 

where ip{oj,rf) gives the probability that a parent with the parameter value r\ produces offspring with 
parameter value uj. Model ()5.5[) . however, is not covered by Theorem [TJ We also note that a very 
similar to (|5.4p model was considered in [22], where, to solve the system, it was conjectured that the 
equation for the mean parameter value is (cf. (|2,13p ) 

m = -/<r 2 (o). 

This means that is was implicitly assumed that the susceptibility distribution is normal, because the 
normal distribution is the one that satisfies the condition o~ 2 (t) = cx 2 (0). 

Finally, all the models we discussed so far are deterministic. There is long history of studying 
stochastic SIR epidemics with multiple classes of susceptible and infectious individuals, e.g., (TJ (H 
[TOT 157] . which is usually centered around the asymptotic distributions of the final epidemic size. We 
announce here that some of the methods presented in this text can be used for studying stochastic 
models [53] . 
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